Advanced Interacting Sequential Monte Carlo 
Sampling for Inverse Scattering 



F Giraud^ 2 P Minvielle^ and P Del Moral^ 

1 CEA-CESTA, 33114 Le Barp, France 

^ INRIA Bordeaux Sud-Ouest, Domaine Universitaire, 351, cours de la Liberation, 
33405 Talence Cedex, France 

E-mail: f rancois . giraudOens-cachan . org 

E-mail: pierre.minvielle@cea.fr 

Abstract. The following electromagnetism (EM) inverse problem is addressed. It 
consists in estimating local radioelectric properties of materials recovering an object 
from global EM scattering measurements, at various incidences and wave frequencies. 
This large scale ill-posed inverse problem is explored by an intensive exploitation of 
an efficient 2D Maxwell solver, distributed on high performance computing machines. 
Applied to a large training data set, a statistical analysis reduces the problem to 
a simpler probabilistic metamodel, on which Bayesian inference can be performed. 
Considering the radioelectric properties as a hidden dynamic stochastic process, that 
evolves in function of the frequency, it is shown how advanced Markov Chain Monte 
Carlo methods, called Sequential Monte Carlo (SMC) or interacting particles, can take 
benefit of the structure and provide local EM property estimates. 



1. Introduction 

Inverse scattering is a topic of major importance; it encompasses various applications 
[H [2] in acoustics, optics and electromagnetism, e.g. medical imaging, tomography, 
ionospheric sounding or SAR (Synthetic Aperture Radar). In electromagnetism (EM), 
the direct scattering problem is the determination of the scattered field, due to the 
scattering of an incident wave in the presence of inhomogeneities, when the geometrical 
and physical properties of the scatterer are known. Conversely, inverse scattering is 
defined as "inferring information on the inhomogeneity from knowledge of the far-field 
pattern..." |2]; it is an inverse problem. In this paper, we focus on a specific, though 
worthwhile, EM inverse scattering issue. The aim is to estimate the electromagnetic 
properties of materials from global microwave scattering measurements. Related 
applications can be located at the crossroads of non-destructive testing, quahty control 
and material measurement. Many EM material characterization techniques have been 
developed in the domain of agricultural and food materials, radar absorbers ^, etc. 
Most of these techniques, from the transmission lines to the admittance tunnel method, 
require small-scale material test samples. For instance, transmission lines enclosed 
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samples inside the conductors of a transmission-line sample holder. Although the 
EM properties (i.e. permeability and permittivity) can be measured, they can differ 
significantly from the final product's ones, when the materials are assembled and placed 
on the full-scaled object or system [3]. The so-called free-space RCS (Radar Cross 
Section: scalar that quantifies reflectivity) methods [3] can overcome this pitfall by 
measuring the monostatic reflectivity of a large planar sample. The sample is then 
located inside an anechoic chamber, in the far field of the transmitting and receiving 
antennas. The reflectivity is measured at various arrival angles of the incident wave. 
Besides, let mention the classic bistatic alternative in near field, known as the NRL 
arch method [5] . In this paper, we focus on the following challenging inverse scattering 
problem: the control and evaluation of EM properties of a full-scaled objet or mock- 
up from the global reflectivity measurements in a free-space RCS device. Deviations 
of microwave properties, such as permeability and permittivity, are to be determined 
along the object. 

Nearly 50 years ago, a closely related issue was formerly outlined in [1]. Least- 
square optimization was applied to determine the dielectric constants that made the 
analytically computed RCS fit with measurements. This issue reemerged in a slightly 
different way in [5j ; both complex permittivity and permeability of a lossless plane 
stratified medium were evaluated. More recently, [6] considers the reconstruction in 
microwave tomography of the dielectric properties of a strongly inhomogeneous object 
by a stochastic global optimization algorithm, based on simulated annealing. Similarly, 
p'J develops a pseudoinversion algorithm for 2D imaging, with the aim locating and 
estimating the dielectric permittivities of unknown inhomogeneous dielectric cylindrical 
objects. On the whole, inverse scattering is known to be an ill-posed inverse problem. 
Like image reconstruction and many other imaging inverse problems [Hlin], it necessitates 
at some step a regularization procedure: it tends to eliminate the artificial oscillations 
resulting from to the problem ill-posedness. According to [2J, the procedures can 
be partitioned into the next two families: the non-linear optimization schemes and 
weak scattering linearization approximation methods, such as physical optics and Born 
approximation. Besides, let mention the efficient linear sampling method in 3D shape 
reconstruction of obstacles due to local inhomogeneities [TOl [2] . 

In this paper, a global statistical approach is developed to address the "free space 
RCS" inverse scattering and solve the large scale ill-posed inverse problem. In some 
way, the approach can be considered to be part of the two aforementioned procedure 
families. It involves an approximation method. Intensive Maxwell solver computations, 
distributed on high performance computing (HPC) machines, results in a surrogate 
likelihood model. It is the starting point of a complete statistical dynamic model 
framework that leads to an efficient inference scheme, close to optimization. It stems 
from statistical signal processing and advanced Monte Carlo sampling (e.g. Markov 
Chain Monte Carlo). Bayesian inference is performed by a sequential Monte Carlo 
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(SMC) stochastic algorithm. These algorithms are called "interacting particles" [TT] 
or particle filtering in adaptive filtering and sequential estimation. They are used to 
provide, in addition to microwave properties estimates of materials, the very significant 
information of the associated uncertainties. From the seminal work of Geman and 
Geman [12], stochastic methods have been commonly used in inverse scattering and, 
more generally, in image inverse problems: simulated annealing for image reconstruction 
|13j . expectation-maximization algorithm for radar imaging [13], etc. In microwave 
imaging, [15] points out genetic algorithms and stochastic heuristics, such as differential 
evolution methods, memetic algorithms, particle swarm optimizations, ant colonies, 
etc. In short, many attempts have been made in electromagnetism to apply stochastic 
methods to tricky inverse problems or non-convex optimization (such as multilayered 
radar absorbing coatings [TU] ITT]). Though powerful, stochastic inverse methods often 
come up against high-dimensional curse. In the approach, it is taken advantage of the 
problem structure to achieve a Rao-Blackwellisation strategy [HI [11] of Monte Carlo 
variance reduction and design a powerful stochastic inversion method that overcomes 
the high-dimensional obstacle. 

This paper is organized as follows. In section|2| free-space RCS material measurements 
is introduced and the inverse scattering problem is developed. Next, section [3] describes 
the probabilistic modeling , from the surrogate likelihood model to the overall statistical 
dynamic modeling framework and, at its core, a hidden Markov model (HMM). 
The inversion Rao-Blackwellised stochastic algorithm is developed in section |4j It is 
evaluated in section |5] where its statistical performance is assessed. 



2. The inverse scattering problem 

2.1. Electromagnetic scattering measurement 

EM scattering measurements have been achieved ever since radar invention [3]. Briefly 
speaking, EM scattering is the standard phenomenon that occurs when an object is 
exposed to an EM wave and disperses incident energy in all directions (scattering is this 
spatial distribution of energy). Some energy is scattered back to the source of the wave. 
It constitutes the radar echo of the object, the intensity of which results from the radar 
cross section (RCS) of the object. More precisely, RCS is defined by: 

as= lim 4rf^5^^^ (1) 
R^+^ lEincI 

It quantifies the scattering power of an object, i.e. the ratio between the scattered power 
density Escat at the receiver and the power density of the incident wave at the target 
(with R the radar-object range). It depends on the wave polarization and frequency. 
The AttR'^ term takes into account the radiated spherical wave. Implicitly, ([T]) requires 
that the incident wave is planar {R — )■ +oo). Practically, it is possible to measure the 
RCS at limited ranges with a sufficient accuracy. It is usually achieved in indoor RCS 
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test chambers, also called anechoic chambers. There, interferences can be limited by 
microwave absorbing materials (see figure [T]). 




Figure 1. RCS measurement inside an anechoic chamber 

In the article, we consider that an object or mock-up is illuminated by a radar, i.e. 
a single antenna or a more complex device (such as the antenna array of figure Q that 
fulfills to a certain extent directivity and far- field conditions [12]. Herein the radar 
system is monostatic, which means that the transmitter and receiver are collocated. 
Its common principle is described in figure [2j Considering that the radar illuminates 
the object at a given incidence with a quasi-planar monochromatic continuous wave 
(CW) of frequency / (incident electric field Einc), the object backscatters a CW to 
the radar (scattered electric field Egcat) at the same frequency. With an appropriate 
instrumentation system (radar, network analyzers, etc.) and a calibration process, it is 
possible to measure the complex scattering coefficient, which can be roughly defined by: 
S = It sums up the EM scattering, indicating the wave change in amplitude and 

phase. S is closely linked to the RCS, with: = It is important to notice that 

the scattering coefficient quantifies a global characteristic of the whole object-EM wave 
interaction in specific conditions (incidence, frequency, etc.). It is possible to measure 
the scattering coefficient for different transmitted and received polarizations. 



rotate 



RADAR 




^ , III III ' ^ 

^scat 



Figure 2. Monostatic scattering measurement principle 
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Let assume the following conventional RCS acquisition mode, widely used in Inverse 
Synthetic Aperture Radar (ISAR) imaging. It consists in measuring various complex 
scattering coefficients S: 

- at different wave frequencies: / G {fi, f2, ■ ■ ■ , fxf}, for Kj successive discrete 
frequencies. Basically, it consists in a series of transmitted narrow-band pulses, 
commonly known as SFCW (Stepped Frequency Continuous Wave) burst pO]. 

- at different incidence angles: 9 G {6^1, ^2, ■ ■ ■ , ^Kq}, for different incidence angles 
(object rotation with a motorized rotating support). 

- at different (transmitted and received) linear polarizations: pol G {HH,VV}, 
meaning respectively, horizontally and vertically polarized both at microwave 
emission and reception. 

Let call Ai the complete measurement, set of 2 ■ Kf ■ Kg elementary complex 
scattering coefficients: M = {>S-^'^'P°'}, for / G {/i, ■ ■ ■ , /x^.} x 6* G {6*1, ■ ■ ■ , x pol G 
{HH, VV}. 

2.2. Nondestructive testing 

In this article, we are interested in an industrial control issue, that can be assimilated 
to nondestructive testing (NDT). Unlike usual EM material characterization techniques 
[5], the point is to determine or check radioelectric properties (i.e. relative dielectric 
permittivity and magnetic permeability) of materials that are assembled and placed on 
the full-scaled object or system. Is it possible from the above complete measurement 
M.1 Is it possible to extract some local information on the material properties along 
the object from the global scattering measurement information? 

area M 




Figure 3. The object coated by Na material areas 

In order to circumscribe the investigation, the article is restricted to a metallic 
axisymmetric object, which is is coated by A^^ material areas, each area corresponding 
to a rather homogeneous material, with its associated isotropic radioelectric properties 
weakly varying within the area. It is illustrated in figure |3| with an ogival shape 
taken from the RCS benchmark [2lj. Consequently, the aim is to determine from 
the global scattering measurement M. the unknown isotropic local EM properties 
(ei, /ii), (e2, /i2)5 ■ ■ ■ ) (e^V! /^Af) along the object, where A^ is the number of different 
elementary zones (cf. Figure |4]). 
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Figure 4. Elementary mesh zones 



2.3. An inverse problem for Maxwell's equations 

Naturally, there is no direct model that is able to compute the radioelectric properties 
from global scattering information. On the contrary, the forward scattering model based 
on the resolution of Maxwell's equations can determine the scattering coefficients given 
the EM properties, the object geometry and acquisition conditions (i.e. wave frequency, 
incidence, etc.). It lies in the resolution of Maxwell's equations, partial derivative 
equations that represent the electromagnetic scattering problem of an inhomogeneous 
obstacle. It is performed by an efficient parallelized harmonic Maxwell solver, an exact 
method that combines a volume finite element method and integral equation technique, 
taking benefit from the axisymmetrical geometry of the shape [22]. Discretization 
is known to lead to problems of very large sizes, especially when the frequency is 
high. Furthermore, as it is shown further, the solver is to be run many times for the 
inversion purpose. Hence, it necessitates high performance computing (HPC): a massive 
supercomputing system, with nearly 20,000 processors and a performance higher than 
1 petaflops (million billion operations per second). 
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Figure 5. The inverse scattering problem 



Figure [5] sums up the entire inverse scattering problem. On one hand, the RCS 
measurement process, that includes acquisition, signal processing, calibration, etc., 
provides the complex scattering measurement Ai, with uncertainties. On the other 
hand, it would be useful to "row upstream" the Maxwell solver, in order to determine 
the unknown radioelectric properties, denoted by x. Yet, even with recourse to HPC, 
there is no direct way to solve what turns out to be a high dimensional ill-posed 
inverse problem, like imaging inverse problems Next, we propose a global statistical 



Advanced Interacting Sequential Monte Carlo Sampling for Inverse Scattering 



7 



inference approach, which is able to take into account prior information and achieve 
the required inversion. Like Tikhonov regularization, it tends to ehminate artificial 
oscillations due to the ill-posedness of the problem. 

3. The statistical problem formulation 

The global statistical approach is introduced gradually, from its formulation at a given 
frequency fk to the whole stochastic model at the various frequencies fi, f2, - ■ ■ , fxf 

3.1. The problem statement at a single frequency fk 

Consider a given frequency fk of the SFCW burst. Let define the two main modeling 
components at fk'- the system state x^, the observation and the probabilistic link 
between them, i.e. the likelihood model p{yk\^k)- To lighten the notations, they are 
denoted respectively x, y and p(y|x) in this section. 

3.1.1. System state x = [e' e" fj,' i-i"]'^ includes the relative permittivity and 
permeability components of the N elementary zones, where ' and " denote respectively 
the real and imaginary parts [f] (at frequency fk)- The four components can be developed 
as: e' = [e'^ ■ ■ ■ e'^f , e" = [e" ■ ■ ■ e'^^f , ^ = [n[ ■ ■ ■ fi']yf and fj/' = [fi" ■ ■ ■ x is in a 
system space of dimension AN; it includes all the unknown parameters that are to be 
estimated. 

3.1.2. Observation y = [3?(iShh) ^^('5hh) '^{<Svv) ^i<Svv)]^ contains the real 
(3^^(•)) and imaginary (Q'(-)) parts of the complex scattering coefficients iShh et iSvv 
measured at the Kq angles 9i,---,9Kg (at frequency fk)- The two complex terms 
Sun and 5vv can be detailed: Sun = [5^^'^^'™ 5^*'^^'™ ■■■ and 
5vv = [5^"'^''^^ S^"'^^'"^^ ■■■ 5/*'^ife'VVj^_ The observation space dimension is 

3.1.3. Likelihood model p(y|x) describes the probabilistic relation between the system 
state X and the observation y (at frequency fk)- In other words, it provides the 
probability distribution of the observation y given a known system state x. It is a 
key element of the knowledge that needs to be taken into account. Our inference goal 
is going to inverse this statistical relation. The likelihood model can be expressed as a 
multidimensional Gaussian of mean Jiv[axwcii(x) and covariance matrix R^: 

y |X ~ A/'(J'Maxwcll(x), Rm) (2) 

where J-iviaxweii is the direct model, from the state space to the observation space, 
that relies on the aforementioned Maxwell solver. Taking into account measurement 
uncertainties, the likelihood model results from the following considerations. 

I In other words, e = e' + je" and /j, = fi' + j/j," (for time dependence convention e-'"'). 
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- The Maxwell solver, based on a direct method, is exact, i.e. extremely precise. 
J-Maxwcii is assumed to compute the "perfect observations", meaning without 
measurement noise, bias, etc. Implicitly, it is assumed that the shape object is 
perfectly known and that, conditionally to radioelectric properties, uncertainty only 
comes from measurement. 

- From previous measurement uncertainty analysisis (see metrology guideline |23j). 
it has been shown that the measurement uncertainty can be reasonably modeled 
by an additive Gaussian noise (y = Jiviaxweii(x) + Vm, ~ A/'(0,Rin)) with the 
quantified covariance matrix Rm- 



Consequently, the likelihood model can be expressed as (with u 
p(y|x' 



g— |(y^-^Maxwcll(x))'^R """(y — -^Maxwell (x)) 



(3) 



(27r)VdetR 

At first sight, just considering a single frequency fk, numerous evaluations of p(y|x), 
i.e. of the Maxwell solver Jiviaxweii(x), are required in order to solve the inverse problem; 
they can be far too time-consuming, even with high performance computing. To avoid 
heavy J-Maxweii computations, a statistical learning approach has been achieved. Its 
basic principle is to build a surrogate model, i.e. an approximation of J-iviaxwcU that is 
acceptable in the limited domain of interest. In a way, it is related to weak scattering 
linearization approximation methods of [2] in inverse scattering, and among them, 
the widely used Born approximation [T, "2]. Here, the statistical linearization is not 
performed from truncation of physical interactions, but from full Maxwell solution 
computations that take into account multiple interactions, creeping waves, etc. The 
system, i.e. the high dimension state space of x and the associated system response 
•^Maxwell (x), is explored by random sampling, according to a prior knowledge about 
the expected radioelectric properties (prior distribution p(x)). The computations are 
massively distributed on HPC machines, each computation involving the parallelized 
Maxwell solver. The computation number depends mainly on the state space dimension. 
The Monte Carlo simulation process leads to the following training set: 

S = {(xW,y«),(x(^),y(^)),...,(x(^^),y(^^))} (4) 



where x'^'^^ ~ p(x) (~ for realization of) and y'^'^^ = J-'Maxwcii(x*^'^-') (for k = 1---Ns), 
Ns being the number of samples. Multidimensional linear regression provides a 
straightforward and efficient way to build a linear model y = /(x) + vi (vi is an 
linearization error term) with: 



/(x) = A.x + y%r/(x) = A*.[l x] , A'^ = [y^ A] 



(5) 



A* is the least square (LS) estimates of the matrix of parameters that minimizes the 
errors to linearity {6i), is given by the solution to the normal equations: 



A* = (A'J ■ XB)~^X^yB with Mb 



1 


x(i) " 




y(l) 


1 


X(2) 


,yB = 


y(2) 








1 









(6) 
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where is the (4A^ x A^^) input matrix and the ( AKg x Ns) response matrix, from 
the training set B. For numerical stabihty, a QR decomposition of is introduced. 
By residual analysis, it is then possible to assess the linear model fitness, i.e. to 
determine the discrepancy between the data and the model in the domain of interest. In 
principle, the covariance matrix (R^) evaluation of the linearization error vi may require 
a supplementary data set or cross-validation methods. Remark that additional statistical 
analysis can be achieved to extract reduced models, removing useless explanatory 
variables, i. e. permittivity or permeability components of zone subsets. That depends 
on the wave interaction, especially on the frequency band. 

Back to the likelihood model (|2]), it leads to an overall error term v = vi + 
of covariance matrix R [§] and to the following linear Gaussian (LG) likelihood model 
(reintroducing the subscript k for frequency fk)'- 

yfc|xfc ~ M{Ak ■ Xfc + y°, R^) or = [A^ ■ + y°] + Vk (7) 

with Afc and learned from the training set Bk- It is illustrated in figure [6] for the 
ogival shape of figure H (A^ = 137, / = 1.5 GHz, ^ = 0° : 1° : 180° - exploration: 
1000 HPC J-Maxwcii simulations). Inside each bloc, the pattern can be explained by the 
coherent contribution of each elementary zone. 




Figure 6. Matrix A^, illustration 

3.1.4- Bayesian approach If such an inversion at a single frequency fk could be solved 
by classical regularization methods [9], Bayesian estimation offers a convenient and 
powerful framework. Let us probabilize the unknown state vector x^. and consider a 
prior probability distribution p(xfc). It is possible to model the priori knowledge with a 
Gaussian distribution: x^ ~ A/" (nifc, P^). 

§ In our context, the linearization error turns out to be much lesser than the RCS measurement 
uncertainties: R + R; ~ R. 



Advanced Interacting Sequential Monte Carlo Sampling for Inverse Scattering 



10 



The mean irik (dimension A^) defines the reference radioelectric properties for the Na 
areas that divide the object (cf. figure [s]). 

-I T 



nifc 



mi 



rrii, 



m7 



m7 



where mf 



4(2) ■■ -61(2) 



(8) 
being 



area i area 2 

the reference real permittivity of area i {i 
mt and m(' . 



area Na 
Similar construction for m' 



k ' 



The covariance Pk (dimension x A^) quantifies the prior uncertainty around nifc. 
Pfe is block-diagonal: = diag(P^', P^", P^ ^Pfc )■ 1^ means that the properties 
(e', e", /i', fi") are assumed to be uncorrelated. Each property block is block- 
structured itself. For instance, P^^ = diag(P^fc (1), (P^^ (2), ■ ■ ■ , (P^ {Na)), expressing 
the assumed property independence between areas. Focusing on one block Pl{i), 
a squared exponential covariance expresses the spatial homogeneity (of the given 
property) between components, i.e. elementary zones of the object that belong to 
the same i*^ material area : 



(Jl 
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■■• Ps 

ps 1 



(9) 



with [cr^'(«)] is the spatial variance of i^^ area and ps € [0,1] the normalized 
spatial correlation parameter (e.g. ps = 0.95). With this Markovian property, 
commonly used in Gaussian field modeling, correlation decreases geometrically with 
the distance between components. Similar construction for P|, , P^ 



P^^ andP^ 



With linear Gaussian structure, i.e. Gaussian prior and linear Gaussian likelihood, 
Bayesian inversion can be performed straightforwardly, with closed-form solutions. In 
our problem, it is a piece of the more complex global problem that encompasses the 
frequency variation. 



3.2. The global problem statement 

Radioelectric properties are known to vary in function of the wave frequency They 
can be quite different from the lower band frequency /i to the higher band one fx- 
The basic idea is to maintain the former statistical modeling at each frequency fk while 
introducing additional a priori information about the dynamic in frequency, i.e. how 
quickly can a property move with frequency, how correlated are a property at two 
different frequencies, etc. This regularity information can be quite different from one 
EM property (e', e", p") to another, as well as from one material to another. 
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3.3. Generalized Auto-Regressive random process 

The statistical modehng extension consists in modehng the whole sequence (x/j,/c G 
{1, . . . , Kf}) by a generalized autoregressive (AR) random process: 



Xl 

Xfc+l 



A/" (mi, Pi) 
rrifc+i + Bp 



H 



k+l 



X,. - m,,) + Jid - D2 . Hfc+i • Vfe (10) 



where is the square root of the covariance matrix P^ [jj] (V^, k E {1, . . . , K}) 
are i.i.d. (independent, identically distributed) M{0,ld) and Dp is a positive diagonal 
matrix commuting with H^. The dynamic model expresses the linear Gaussian 
correlation structure. It can be checked that the marginal distribution of x^ is still 
More generally, it can be shown that the distribution of concatenated 
,^Kf) is Gaussian with mean m = (mi, . . . , m/^',) and covariance 



vector X = 
matrix: 



Xl, 



Id 

Dp 
D^ 



-L'p 



Dp 
Id 
D„ 



D^ 
D„ 



-L'p 



D, 



D. 



'111 



where H is the block diagonal matrix 7i 
joint distribution (xj,Xj) is expressed . 



Id 

diag(Hi, 



,Hj^^.). Basically, every 



The matrix Dp takes into account the frequential correlations of the EM properties 
Xl ■ ■ ■ x^^ ; it refers to a hyper-parameter p. According to about frequency correlation 
prior knowledge, the following alternatives can be considered: 

(i) The frequency correlation doesn't depend on the material and the EM property (e', 
e", jj! or ji"): p is scalar (g [0, 1]) and Dp = p.Id- 

(ii) It depends on the material: p is A^^a-dimensional (g [0, 1]^"), and Dp is the block- 
diagonal matrix made up of Na terms pj.Id- 

(iii) It depends on both: p is 4.A^a-dimensional and Dp is the block-diagonal matrix 
made up of i.Na terms Pi-Id- 

3.4. A conditionally hidden dynamic Markov process 

The generalized AR random processes include the linear Gaussian models at the various 
frequencies fk {k = 1 ■ ■ ■ Kf). It provides a spatial and frequential correlation structure. 
Assuming that the material areas are known to be quite homogeneous, the spatial 
correlation parameter can be fixed (typically ps = 0.95). Quite the reverse, frequency 

II unique symmetric definite positive matrix such as: • = P^. 
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correlations can not be really known; they are to be determined by the inversion process. 
Back to Bayesian statistics, it is chosen to probabilize the unknown hyper-parameter p. 
Finally, the combination of the AR dynamic model (11) and the likelihood model ([T]) 
end in the following state-space model, observed at "times" fk {k = 1, ■ ■ ■ , K): 

Xfc+i = ■ Xfe + Wfc Yfc = [Afc ■ Xfe + y^] + Vfc (12) 

assuming the initial state xi ~ A/'(mi, Pi). is a transition matrix and a Gaussian 
model noise (E(wfc) ^ 0). Both directly arise from (11); they are not detailed here for 
clearness. 




Figure 7. A graphical representation 



Again, let emphasize that the dynamic model involves that each marginal complies 
with Xfc ~ A^(mk, Pk). On the other hand, it is important to remark that, conditionally 
to the frequential correlation parameter p, the model is a classic linear Gaussian hidden 
dynamic Markov process. A graphical representation of the entire model is given in 
figure [7} Given a value of p, the lower part describes a linear gaussian system. The idea 
is to make the most of this specific structure. 



4. Advanced Sequential Monte Carlo inversion 
4.I. The Rao-Blackwellized Approach 

As already mentioned, the unknown hyper-parameter p is probabilized, so it is given a 
prior distribution p{p), assumed calculable (up to a normalizing constant) and easy to 
sample. The posterior distribution p(x, p|y) can be decomposed as: 

P(x,p|y) =p(x|p,y) ■p(p|y) (13) 

Since the system is linear Gaussian conditionally to p, the conditional distributions 
p(xfc|p, y) can be straightforwardly computed by classic Kalman filtering. This forward 
algorithm can be completed by backward smoothing, in this off-line context; the overall 
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is often called "Kalman smoother". On the other hand, the term p(p|y) can been 
decomposed as: 

P(p|y) oc p{p)-p{y\p) 

Kf 

(14) 



p{p)-T[p{yk\p,yi,---,yk-i] 



^ ^ ■=Jk(p) 



Again, for any hyper-parameter p, the quantities Jfc(p) can be evaluated from the 
likelihood terms provided by the Kalman filter. Eventually, it is possible to exploit 
this conditional system structure, with Kalman smoothers that can be applied and 
integrated in the following interacting particle approach. In a first step, a stochastic 



algorithm (described in section 4.2) gives an approximation of p(p|y). It estimates the 
frequential correlations (i.e. regularity) of the EM properties e'(/), e"(/), p'(/), p"(/). 
In a second step, the first moments of can be evaluated (for each frequency fk) by 
the theoretical conditioning relations: 

E(xfc|y)=E[E(xfc|p,y)|y] (15) 

Var(xfcly) = E[Var(xfe|p,y)|y]+Var[E(xfe|p,y)|y] (16) 

Note that Kalman recursions are used both in the first step for calculating the likeli- 
hood of the hyper-parameter p (up to a normalizing constant) and in the second step 
for determining the quantities E(xfc|p, y) and Var(xfc|p, y). This idea of mixing analytic 
integration (here Kalman evaluation of p(x|p, y)) with stochastic sampling (here to ap- 
proximate p(p|y)) is a variance reduction approach, known as Rao-Blackwellisation 



Let us denote by ri{dp) the probability measure associated with the marginal 
distribution p(p|y), for a fixed observation vector y. Similarly to [18J, we choose to 
implement for the first step an efficient interacting particle approach, called Sequential 
Monte Carlo (SMC), in order to estimate rj. We now give a brief but general description 
of these methods. 



4.2. The SMC algorithm 

Sequential Monte Carlo is a stochastic algorithm to sample from complex high- 
dimensional probability distributions. The principle (see, e.g., [TT]) is to approximate 
a sequence of target probability distributions (r7„) by a large cloud of random samples 
termed particles {Cn)i<k<Np £ E^^ , E being called the state space. Between "times" 
n — 1 and n, the particles evolve in the state space E according to two steps (see figure 
§: 

(i) A selection step: every particle Cn~i is given a weight Ui defined by a selection 
function Gn '■ E ^ (0, +oo) (i.e. Ui = Gn{Cn-i))- resampling (stochastic or 
deterministic), low- weighted particles vanish and are replaced by replicas of high- 
weighted ones. 
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(ii) A mutation step: each selected particle Cn-i moves, independently from the 
others, according to a Markov kernel M„ : E ^ E. 



Gn 



selection 



S-1 



5n— 1 



Mn 



Sn 



mutation 



Figure 8. The SMC 2-step evolution 



Evolving this way, the cloud of particles, and more precisely the occupation distribution 
Vn'' '■= J2k=i "^Cfe (sum of Dirac distributions), approximates for each n the theoretical 
distribution ry„ defined recursively by the Feynman-Kac formulae. It is associated with 
the potentials Gn and kernels M„ (see [21] for further details). More precisely, this 
sequence //„ is defined by an initial probability measure r/o and the recursion: 

r]n = '^G„{Vn-l).Mn (17) 

where "^CniVn-i) is the probability measure defined by '^G„{Vn-i){dx) oc Gn{x).rin-i{dx) 
and, for any probability measure /x, /i.M„ is the measure so that /x.M„(A) = 
J^Mn{x,A)fi{dx). 

The SMC approach is often used for solving sequential problems, such as filtering 
(e.g., flEl |2ni [2Z])- In other problems, like ours, this algorithm also turns out to be 
efficient to sample from a single target measure 77. In this context, the central idea 
is to find a judicious interpolating sequence of probability measures (r7„)o<fc<n/ with 
increasing sampling complexity, starting from some initial distribution 770, up to the 
final target one r/^^ = rj. Consecutive measures rjn and rjn+i are to be sufficiently 
similar to allow for efficient importance sampling and/or acceptance-rejection sampling. 
The sequential aspect of the approach is then an " artificial way" to solve gradually the 
samphng difficulty. More generally, a crucial point is that large population sizes allow to 
cover several modes simultaneously. This is an advantage compared to standard MCMC 
(Monte Carlo Markov Chain) methods that are more likely to be trapped in local modes. 
These sequential samplers have been used with success in several application domains, 
including rare events simulation (see |2H]), stochastic optimization and, more generally, 
Boltzmann-Gibbs measures sampling ([29j). 



4-3. Interpolating sequences of measures 

Back to our objective of sampling from ri{dp), let us denote by E the state space of 
the variable p (i.e. E = [0,1], [0,1]^" or [0,1]"^^"). We have to define a sequence of 
distributions {rin)o<k<nf from the initial distribution rjoi^dp) = p{p)dp (easy to sample) 
to the target one rjnAdp) = rj{dp) = p{p\y)dp. 
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4-3.1. The guiding principle With this in mind, we first define an interesting class of 
Markov kernels on E: let h he a. positive, bounded function on E, and let Q{x, dy) 
be a Markov kernel on E, assumed reversible w.r.t. the Lebesgue measure on E. The 
Metropolis-Hastings kernel Kh,Q{x, dy) associated with h and Q is given by the following 
formula: 

Kh,Q{x,dy) = Q{x,dy). mm (l,'-^^ Wy x 

Kh,Q{x, {x}) = 1 - y <5(^' dy). mm (^1, ^ j 

Using an acceptance/rejection method, this kernel is easy to sample as soon as one can 
sample Q{x,dy) and calculate the ratios h{y)/h{x). Here is a crucial property: if 
denotes the probability measure defined by f^h{dp) oc h{p)dp, then it is well known (see, 
e.g., [30]) that Kh,Q admits ph as an invariant measure: 

More generally, this property is satisfied for the iterated kernel KJ^q, i.e. ph-KJ^Q = ph 
(for any integer m). 

Let rjn be a sequence of probability measures defined with some positive, bounded 
functions /i„ so that: T]n{dp) oc hn{p).dp. Then, for any sequence of reversible Markov 



kernels Qn and any sequence of integers rrin, f]n satisfies the Feynman-Kac formula (17) 
with potentials (?„ := hn/hn-i and Markov kernels Mn ■= KJ^^q^ {Kh„,Q„ iterated m„ 
times). Practically, the consequence is that such a sequence ?7„ can be approximated 
using a SMC algorithm as soon as one can calculate the functions hn up to a normal- 
izing constant. Similarly to traditional MCMC or simulated annealing methods, this 
algorithm is all the more robust when the iteration numbers m„ are large, since the 
kernels Kh^ q^ are just defined and used to stabilize the system. 



4.3.2. Design of bridging measure sequences From these considerations, we propose 
three scheme variants of interpolating sequences of measures. 

(i) The annealed scheme: the sequence rjn is defined by the positive, bounded functions 

hn{p) = Piylp)"" ■ p{p) 

where (an)i<n<n/ is a sequence of numbers increasing from to 1 (arbitrarily 
chosen). In this situation, the potentials Gn{p) used in the selection are equal to 
p(y |p)""~°"~^- Thus, a„ is to be chosen to control the selectivity of these functions, 
which is important in practice. Annealing or tempering is frequently used in SMC 
(see [311 [11]); it is related to simulated annealing (with inhomogeneous sequence of 
MCMC kernels). 
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The data tempered scheme: for all n E {0,1, . . . , Kf}, rjn is the probability measure 

n 

associated with: hnip) = p{p) ■ |Tp(yfe|p, Yi, • • • ,yfc-i)- In other words, at each 



=1 ^ 



generation n, the selection potential Gn{p) that is applied to the particles is the 
term p(y„|p,yi, • . . ,y„_i), i.e. the likelihood of the n-th observation vector given 
the previous ones. This allows the algorithm to work "online", since it treats 
sequentially the observations. According to [3T], it is efficient for problems that 
exhibit a natural order (e.g. hidden Markov models). Yet, when these potentials 
turn out to be too selective, the SMC algorithm turns out to perform poorly since 
the cloud of particles loses its diversity at each selection step. It is substituted for 
the next scheme that overcomes this drawback, 
(iii) The hybrid scheme: similarly to the previous one, this scheme incorporates the 
observations one after the other, but each likelihood function Jfc(p) is handled as a 
product: 

n, (k) (k) s 
MpY'^' 

i=l 

where for all k G {l,...,Kf}, (ai^^)i<i<nfc is a sequence 0/^1. Then, if 
n = (rii + ■ ■ ■ + Ur-i) + s, the function hn is given by: 



hn{p)=p{p)-[l[jk{p)]-Jr{p) 



(r) 



\k=l 

Note that the selection potential Gn = Jr " " can be arbitrarily controlled. 

For each of these interpolating schemes, the functions hn are calculable up to a 
normalizing constant (Kalman equations), so that the Metropolis-Hastings kernels 
(possibly iterated) can be used to perform the mutation steps. 

4.4- The global estimation 

To sum up, the joint distribution p(x, ply) can decomposed and evaluated as follows: 

KF output prior 

^p{y\p) -pip) 
p(x,p|y)= p(x|p,y)_ ^ p(pjy) 

KF (+ smoothing) SMC 



As previously mentioned, the SMC algorithm of section 4.2 provides in the first stage 



an evaluation of the frequency correlations p{p\y) (i.e. an approximation fj = rjuf of rj). 
It is computed from the last generation of particles {p^^\ . . . , p'-^''^) := (d^, . . . , C»v )■ 
In the second stage, estimators of EM properties are straightforwardly computed from 
conditioning relations (15) and (16) (see details in annex]?]); it consists in approximations 



of the mean and covariance matrix of the system state x^. Focusing on a given frequency 
or on a fixed zone, the SMC method provides useful information: 
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For any frequency fk, it computes an approximation of the mean and covariance 
matrix of the system state x^. Roughly speaking, one can sample from the posterior 
distribution p('Kk\y) by picking a p^*^ from the final cloud of particles and computing 
associated samples of x^ by a Kalman smoother conditionally to p^*-* (see further 
page 



illustration figure 12 



20). 



For any fixed zone, the method provides estimators of the mean and marginal 
variance for every frequency, so that the results can be presented as frequential 
profiles, with marginal uncertainties (using the diagonal values of E^) (see further 
illustration figure 13 page 21). 



5. Applications 



In this section, the inverse scattering approach is applied to EM scattering measurements 
of a metallic ogival shape object. The validation is achieved with simulated data in a 
wide frequency band from / = 200 MHz to 8 GHz. Section describes the reference 
nondestructive testing scenario. Next, section |5.2| describes the inversion process and 



illustrates some results. A detailed performance analysis is developed in Section 5.3 
Then, in Section 5.4, we briefly analyze some variants of the approach. 



5.1. Nondestructive testing scenario 

The metallic object We consider the metallic axisymmetric object, previously shown 
in figure [3] its ogival shape, derived from the RCS benchmark |2T], is perfectly known. 
The 2 m long object is coated by A^ = 5 material areas, the isotropic radioelectric 
properties weakly varying within each area. For each material area, the true EM 
properties x^j.^g(/) undergo the following model: ^imQif) = ^xefif) + c ■ A(/). 



1.5 
1 









— V 

A , 
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1 1 1 
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Frequency (GHz) 



Figure 9. The functions A 



At each frequency /, the true (unknown) vector x^j^^g(/) is 4 A = 76- dimensional, 
where: 

- x^gf(/) is a reference frequency profile, depending on the area and on the 
radioelectric component (e', e", p', /i"). Note that these 4Aa = 20 reference profiles 
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are chosen regular and with typical orders of magnitude (i.e. non-negative and 
< 20). 

- A(/) is a perturbation function depending on the radioelectric component. Thus, 
the 4 functions A^', A^", A^/, A^// define the perturbation shapes . As shown in figure 
[oj they are chosen more or less regular (in order to test the inversion capabilities). 

- c is a simple scaling factor, depending on the area. To examine the perturbation 
amplitude influence, increasing values of c are chosen: {0.5, 1, 2, 4, 8}, related to the 
5 successiveareas. 



(Simulated) scattering measurements According to the conventional RCS acquisition 
mode described in section 2.1, complex scattering coefficients are measured for 
both polarizations HH and VV, at Kf = 20 regularly spaced frequencies (/i = 
0.2 GHz, ■ ■ ■ , fxf = 8 GHz) and at Kg = 23 regularly spaced incidence angles 
(^^i = 0°,---,^^, = 180°). 

The observation data y = (yi, . . . , yxf) is simulated from the likelihood model ([2|. 
That involves to run the parallehzed harmonic Maxwell solver (J-Maxweii) and to draw 
an additive white Gaussian noise of marginal standard deviation Un = 10"^ (~ 1%). 
Note that each of the 20 observation vectors y^ is 4 x Kg = 92-dimensional. The data 



is represented in figure [TOj Note on the amphtude representations the high specular 
reflections when the ogival object is turned perpendicularly to the wave propagation 
direction. 



Amplitude (dB) 
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Figure 10. Observation hologram, amplitude and phase (polar HH and VV) 
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5.2. Inversion process 

The goal is to estimate the radioelectric properties, the x^^j^^g term function of the 
frequency /, from the scattering measurements. In this section, we give a few 
implementation details regarding the application context. 

State space The state space dimension stems from the wave frequency number and 
from the discretization of the object in elementary mesh zones. In order to limit it, the 
cutting up of the object is here restricted to = 19 elementary zones. 



Prior information The prior information (see section |3]) needs to be detailed in this 
context. Concerning the prior spatial information p(xfc), its means nik are given, for each 
k, by the former reference frequency profiles x^g£(/fc). Around them, the uncertainties 
are given by the block-structured covariance matrices Pfc of ^ with: ps = 0.95 and 
crfc(z) = 1+0.15 X mfc(i) for any elementary zone i. In other words, we assume a minimum 
standard deviation of 1 that increases proportionally to the reference amplitude value. 
Regarding the prior frequential information, we assume that p depends on both area 
and EM property (e', e", p', p"), so that it is 20-dimensional. As for its prior distribution 
p{p), we set: 



20 



p{p) = n^*^^*) 

i=l 

where all the marginal prior distributions p{pi) are identical and presented on figure 



11 Note that this distribution p(p) can be sampled straightforwardly by sampling 



independently each component pi using, e.g., an acceptance/rejection method. 



Q. 
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Pi 

Figure 11. Marginal prior distribution p{pi) 



Likelihood model The surrogate likelihood model ([7]) has been formerly learned: 
and y° are known (see figure |6]), as well as the marginal standard deviation cxn which is 
in conformity with the measurement noise of the above observation simulation. 



SMC tuning The sequence of probability measures ?7„ is standardly defined by the 
annealed scheme (see section 4.3). To ensure a stable behavior of the SMC algorithm 
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(i.e. keep a good approximation rjn'' — rjn until the end), we chose the following efficient 
adaptive strategies (that make it possible to limit the number of particles to Np = 100): 

- selection step: as mentioned, the increment Aa„ = a„ — q;„_i controls the selectivity 
degree. If Aq;„ is too small, every particle is given approximately the same weight, 
and there is no selection among them. If Aa;„ is too large, the majority of the 
particles are killed, the cloud loses all its diversity, and the SMC algorithm performs 
poorly. Therefore, instead of choosing beforehand Aa„, it is defined adaptively so 
that the selection step kills around 25% of the particle population. This is a way 
to ensure a reasonable selection. 

- mutation step: the mutation step is crucial since it allows the particles to explore 
the state space E. We use Markov kernels M„ defined as being the composition of 
several Metropolis-Hastings kernels Kn^ whose proposition kernels Qn\x,dy) are 
uniform, centered in x, and associated with a window size o"prop,n- To be sure that 
the particles move in a well-sized neighborhood, (i.e. large enough to explore E and 
small enough to converge), the sequence (crprop,n)j always starts with large values 
and decreases geometrically. Once more, we use an adaptive criteria to stop the 
process. 



Results In the context of this reference study, the inversion process takes about 30 
minutes with a current standard processor. Note that the higher the dimension space 
is, the longer the inversion. In figure 12, we show the estimations of /i' for all the zones 
of the object, with their associated uncertainties, compared with the true values, at a 
fixed frequency /14 = 5.6 GHz. Note that the EM property deviation is important in 
our example (see figure |9]). As already mentioned, it is possible to provide some samples 
of the posterior distribution p(xi4|y) to determine the uncertainty on the estimators. 
The EM radioelectric properties are correctly inferred all along the ogival object and its 
5 material areas. The uncertainty recovers more or less the real profiles. 
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Figure 12. EM estimated properties at frequence / — 5.6 GHz 



Figure 13 presents frequential profiles for a fixed elementary zone (the 18* ). All the 



components (e', e", /x', /i") are represented. Each of them is quite accurately estimated. 
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The results are good, even when the perturbations (i.e. the difference between the 
prior and real profiles) are large and irregular. This robustness is due to the adaptive 
estimation of p's components. Next it is confirmed by several thorough analysis. 
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Figure 13. Estimated EM properties of the 18* elementary zone 



5.3. Performance analysis 

To extend the results, we propose a statistical performance analysis of the inversion 



process. It is lead in the same context of section 5.1 As the developed interacting 



particle approach is partly stochastic, two different aspects must be studied. Firstly, for 
a single given data y, the variance of our estimators and Sfc, only due to the random 
feature of the method. Secondly, the average variance of our method for several data 



5.3.1. Stochastic variation For a given data y, our method mainly provides 2 sequences 
of estimators. The posterior mean estimators (xi, . . . , xk/)? and the posterior covariance 
matrices estimators (Si, . . . , Sa-^). As with all stochastic algorithms, one has to verify 
that despite random, it always gives the same result, or at least that its own variance 
is negligible. 

Let X denote the concatenation of the vectors xi, . . . , Xi^^.. Let a denote the 
concatenation of the estimated marginal uncertainties (square root of the S^'s diagonal 
values). Defined in this way, x and a can be considered as 2 matrices of size 76 x 20, 
and the 2 main estimators of our method. 

To quantify the stochastic variance, we simulate an observation data y, and we 
perform the inversion method 30 times. At the end, we get 30 pairs of estimators 
{(x«,a« ),..., (x(30),a(30))}. For any pair of index (i, k) G {1, . . . , 76} x {1, . . . , 20}, 
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we consider the mean values of the estimators and their RMS (root mean square) values: 

^30 -.30 



xu, 



r=l 



r=l 



RMS fx 



RMS (a 



^) ihk) := (^^^(xM(^,A:)-X(^,A:))'j 



1/2 



1/2 



The numerical results, taken over all the pairs of index {i,k), are summed up in 
table [TJ Two points can be clearly emphasized. First, the standard deviation of the 
x*^^) is very small in an absolute way (~ 10"^). Moreover, it is negligible compared 
to the estimated variance of our estimators (at least 1 decade). Secondly, the standard 
deviation of the a^^^ is even smaller (^ lO^'^) and negligible compared to the values of the 
o"*^^^ themselves (at least 2 decades). Consequently, there exists a stochastic variance, 
but it is far negligible compared to the uncertainty inherent to the inverse problem, 
including measurements. 



mean RMS (x) 


max RMS (x) 


RMS(x) 

mean — 


RMS(x) 

max — 


4.16 10-3 


4.11 10-2 


1.08 10-2 


9.87 10-2 


mean RMS (a) 


max RMS (a) 


RMS(d-) 

mean — 


max — 


1.10 10-3 


7.56 10-3 


2.99 10-3 


1.84 10-2 



Table 1. RMS results of x and a 



5.3.2. Average precision The average precision is analyzed on several cases. For this 
purpose, 30 independent observation data {y*-^-',. . . ,y''^''^} are simulated. For each of 
these observation vectors y^^-*, the inversion algorithm computes the pair of estimators 
(x^''\ (T*-^^). The comparison with the true values of x is quantified by the following root 
mean square error (RMSE) : 

RMSE(2,A;) := ( V (xW(2, A;) - Xt,uc(^, A;))' ' 



These made errors are shown on figure [I4j, where they can be compared to the 
estimated errors a. From these results, these conclusions can be drawn. Despite the 
large amplitude and irregularity of the perturbations, far from the assumed prior model, 
the estimators x*^'"-' give a good approximation of xtme (note that the mean RMSE = 
3.68 10-1). Moreover, the RMSE values are comparable to the marginal uncertainties 
given by a, which proves that the estimated posterior variances make sense. 
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Figure 14. RMSE (left) and estimated marginal uncertainties (right) 



In this inversion process, the role of p's estimation is very interesting. Roughly 
speaking, it is as if it can give in advance the shape type of each of the unknown true 
frequencial profile, by estimating its regularity. On figure 15 , we show the results given 
by (x'-^\ 5"'^^^) for the zones number 2, 9 and 17 and the permeability fi'. On the right 
part, the histograms represent the posterior distribution of p. 

As predictable, the difficulty is increasing from zone 2 to zone 17. It is due to the 
perturbation which is larger and larger, as well as irregular. On the right side of 
the figure, we show the histograms of all the particles . . . , p'^^°°)) (each particle 
being represented by its associated component). We clearly see that the more irregular 
is the true signal, the smaller are the p^^\ which is quite coherent since p quantifies 
frequential correlation. Meanwhile, we verify in the center of the figure that, in spite of 
the increasing difficulty, the mean RMSE remains stable. Again, let us stress that the 
adaptive behavior of p estimation is essential to the algorithm robustness. 



5.4- Additional analyses 

We propose now to briefiy analyze the infiuence of other parameters, that can come 
from the context or from the inversion process itself. 



5.4-1 ■ Influence of the Processing Parameters The inversion process we described in 
section |4] admits several qualitative and quantitative degrees of freedom, in the SMC 
step particularly. We propose here our empirical remarks about some of them. 



The number of particles Np Like a classic i.i.d. (independent and identically 
distributed) sampling method, the SMC algorithm precision is proportional to Np^^'^. 
However, in our problem, the main objective is not to have a precise estimation of rj, but 
of X. As the impact of a local variance of p on x is rather small, the crucial point is that 
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Figure 15. n' estimators for zones 2, 9 and 17 



the global cloud of particles reaches the correct area in E. From this point of view, the 
important condition is the stability of the Feynman-Kac flow (see [21]), which ensures 
that the particles don't get lost in E. This is precisely the purpose of the adaptive 
strategies inside the selection and mutation steps). That's why it seems useless (and 
time consuming) to use a high number of particles. Note that below Np ~ 40, the SMC 
approach may be trapped by some local modes. 



The interpolating scheme rjn In addition to the annealed probability measure scheme, 
the hybrid one has been tested. It can assimilate the observations one by one, and 
update the estimators progressively. Moreover, it manages the computational problems 
of selectivity that affects the data tempered scheme. The results are good, nearly 
identical to those obtained with the annealed scheme. And yet, the SMC algorithm 
lasts around 4 times longer than before. 




Generation number n 



Figure 16. Annealing parameter al 



(fc) 



This behavior can be easily interpreted by figure 16 It appears that many 



observations do not bring any new information, so that the associated annealing 

(k) 

sequence a) takes value 1 at once. On the contrary, when a new observation provides 
information in contradiction with the previous ones, the particles have to migrate from 
an area of E to another, which takes a longer time (more steps). 
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The parameter p The prior distribution p{p) of figure 11 has a limited impact on the 
final estimation of rj. Corresponding to a prior knowledge of frequency regularity, it is 
arbitrary chosen in order to penalize the small values and favor regular profiles. But in 
practice, this penalization term p{p) is less determining than the likelihood one p{y\p)- 
Besides, p can be defined 5-dimensional. In this case, the SMC algorithm performs 
quicker. However, the underlying hypothesis, i.e. the frequential correlation is the same 
for e', e", p\ p", is not necessarily fulfilled in practice. 



5.4-2. Context influence As we mentioned, the method is very robust concerning the 
amplitude and the irregularity of the perturbation (deviation from the reference profiles). 

Measurement noise However, it is naturally sensitive to the observation noise 
magnitude. Its performance degrades when the observation noise is too high. That 
is clearly a matter of information. Numerically, it can be explained by considering the 
accurate approximation given by the surrogate model. Indeed, the matrices are 
ill-conditioned. In particular, p' and e" components are highly correlated; it is the same 
for p" and e'. 




Figure 17. Estimation of /i' and e" , an = 10 



In figure 17, we give the estimations of these 2 quantities in the case where the am- 
plitude noise = 10~^. One can then see that the unknown perturbations of /z^rue ^^^1 
e"j.^g are correctly detected by the process, but improperly distributed between p' and e". 



Problem dimension Concerning the computation time, it is very sensitive to the 
dimension of the problem. The reason is simple: each elementary evaluation involves 
(among others) a Kalman smoother, i.e. Kf inversions of 4N x 4A^-sized matrices. 
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Consequently, for a problem of large dimension, it could be appropriate to parallelize 
the SMC algorithm and distribute the Kalman smoothers on HPC. 



6. Conclusion 



An efficient statistical inference approach has been applied. From global EM scattering 
measurements, it manages to estimate local radioelectric properties of materials 
assembled and placed on the full-scaled object. The inverse problem is solved by 
combining intensive computations with high performance computing (HPC), surrogate 
modehng and advanced sequential Monte Carlo techniques dedicated to frequency 
dynamic estimation. It takes advantage of the problem structure to achieve a Rao- 
Blackwellisation strategy of Monte Carlo variance reduction. On top of that, the 
Bayesian approach quantifies the uncertainties around the estimates. 

To tackle higher dimensional problems, it could be interesting to apply close 
stochastic techniques, such as "interacting Kalman filters", and above all, benefit from 
the highly parallelization/distribution potential on HPC to tackle 3D geometries and 
high-dimensional problems. 



7. Annexe: estimation of and conditioning 

For a given y, it is possible to define judicious estimators of x^ (for each k). Indeed, the 
first moments x^ and S^) can be determined from the theoretical conditional expecta- 
tion Xfc := E[xyt|y] and covariance matrix := Var[xyt|y]. 

For all p E E, let set: Xfc(p) := E[xfc|p,y] and Sfe(p) := Var[xfc|p,y], i.e. the 
main quantities provided by the Kalman smoother. Under this notation, we combine 



f} ~ ?7 together with the equation (15), and derive a natural choice for the estimator x^: 



Xfc = E[xfc|y] = E[E(xfc|p,y)|y] 



Np 



1 



1=1 



Regarding the covariance estimator 1]^, under the same notation and according to (16) 
we have: 

Sfc = E fSfe(p)ly) + Var (xfc(p)|y) . 



k 



We estimate S^,^-* and S^^-* separately: 
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(i) Evaluation of S[.^^ 



peE 



^ i=l 



(ii) Evaluation of Tr^ 



(2) 



.(2) 



E 



Xfc(p) - Xfc) (xfc(p) - Xfc) |y 
j {^k{p) - Xfc) (xfe(p) - Xfc)'^ ri{dp) ~ J (xfe(p) - Xfc) (xfe(p) - Xfe)'^ ?7(cip) 



peE 



p i=l 



Finally, the estimator of Sfc is given by: Sfc := + S^^-*. 
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